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A new set of elements is proposed to describe Keplerian motions subject to perturbing forces. The resulting equations 
do not break down for small eccentricities, small inclinations and rotations of the ideal frame reference that are half turns. 
The parameters are selected with a view of simplifying the programming of the right-band members. 
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1 . Introduction 

With the development of minicomputers and micropr cessors, the advantages of the method 
of variation of parameters (in short, VOP, also referred to as the variation of elements or the varia- 
tion of constants), may outweigh, in many individual problems, the simplicity and generality of a 
straightforward numerical integration in Cartesian coordinates with respect to a fixed frame of 
reference (CowelPs method). These advantages include especially (1) a lesser accumulation of 
errors because the integration processes parameters that are affected solely by the usually small 
perturbations, (2) fewer steps of integration because the parameters are slowly varying, (3) a deeper 
insight in the physics of the problem because the elements are selected to enhance the special 
effects of the perturbing forces on the characteristics of a Keplerian motion. 

There have been numerous attempts at finding the parameters best suited to the differential 
equations 

x = X, (1) 

*— (g)* + F (2) 

that describe the motion of a mass particle with position x and velocity X in a Newtonian potential 
of constant fi to which is superimposed a perturbing force F. Classical astronomers have been overly 
concerned that conversion from Cartesian coordinates to Keplerian elements not raise the order of 
the original system. For someone integrating numerically a differential system by consulting a table 
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of logarithms or operating a desk calculator, higher order necessarily meant loss of efficiency. 
The breakthrough in adapting VOP to modern computing conditions was made by Samuel Herrick 
[l] 1 . Abandoning the requirement that the order of the system be kept at its minimum, be developed 
into a computer algorithm the suggestion by Milankovic [2] that the Eulerian angles h, g, and / 
(respectively the right ascension of the node, the argument of the perigee and the inclination) 
be incorporated in the direction cosines of the major and minor axes of the reference ellipse. 
These unit vectors, henceforth denoted a and b, when taken as parameters, produce an overabund- 
ant system, yet they render the right-hand members of the equations more symmetric, hence 
simpler to program and faster to evaluate. Applied to the orbit of Icarus, the new version of VOP 
amply justified Herrick's expectations [3]. 

Vectorial elements had been used by Stromgren [4], Bilimovic [5, 6], Musen [7], and Popovic 
[8] to derive equations for combinations of the classical Keplerian elements. Now following the 
lead taken by Herrick, Musen [9] proposed another computer algorithm in which the direction 
b of the minor axis and the normal to the orbital plane (henceforth denoted h) are used as elements. 
Musen's algorithm has been the root of the computer programs developed at Cincinnati Observa- 
tory in relation with the Vanguard Project and the Minor Planets Center. 

Substituting the couple of directions (a, b) or (b, h) for the Eulerian angles (h, I, g) amounts to 
referring perturbed Keplerian motions to a slowly moving frame. This viewpoint is not new; Hansen 
[10] designed his lunar theory around the concept of an ideal frame chosen to ensure a measure of 
independence to the motion in a plane from the rotation of the intermediate referential. Herget 
[11] seized on this idea to eliminate the singularity caused by small eccentricity in the FOP equa- 
tions based on the apsidal frame (d, b, h). But as he chose to represent the motion of the ideal frame 
as a finite rotation in order to keep the order of the system at its minimum, Herget introduces 
another kind of singularity, namely, those due to rotations that are half turns. We show here how 
Eulerian parameters eliminate the singularities due to small inclinations without generating 
spurious singularities at half turns. 

Wider acceptance of the method of VOP is hampered by the lack of a systematic derivation 
of the various concepts and equations free from second-hand references to foreign journals and 
monographs out of circulation and in a mathematical language familiar to contemporary engineers. 
We have endeavored to restructure the piecemeal contributions brought to the theory of special 
perturbations over the last forty years. In the course of this work, we have come across a new set 
of elements that we think are particularly well suited to numerical integration by electronic com- 
puters. 

2. Ideal Frames and Departure Points 

The angular momentum (per unit of mass) of a particle is the vector 

G = xXX; (3) 

let G and h be respectively its norm and its direction. Let x be the unit vector in the direction of 
&, and y the unit vector such that the triple (x, y, h) makes an orthonormal base. We call it the 
orbital frame. By definition 

G=Gh, (4) 

x= r i, (5) 

f=hXx. (6) 



1 Figures in brackets indicate the literature references at the end of this paper. 



The base (i , yo, n ) that is the orbital frame at the epoch to is adopted as the fixed frame of 
reference. 

In the instantaneous orbital frame, the perturbing force is decomposed in the sum 

F=Rx+Ty + Nh; (7) 

the components R, T and N are sometimes referred to as the radial, the transversal and the normal 
components of the perturbation respectively. 
Differentiation of (3) yields readily that 

G=xXF. (8) 

Then by differentiating the relation 

C* = GG, 

one obtains that 

G=rT, (9) 



d h = ~Ny. (10) 



and on considering (8) and (9) while differentiating (4), one gets finally that 

dt n G 

Likewise, from the identity, 

r 2 = x • x, 
one deduces that 

r=x-X (11) 

and comes, through (1), (2) and (11), to the basic equation 

d . G . 

Finally differentiation of (6) via (10) and (12) produces the last equation 



(12) 



**—+* + G N * ,13) 

necessary to describe the motion of the orbital frame. 

Accordingly, the motion of the orbital frame with respect to the fixed frame may be viewed 
as the rotation 

d m m d . A d „ 

— x=a>Xx,— - y = a) X y, — - n=(oX n 

dt dt J J dt 



with the angular velocity 



io = -n + -N x. (14) 



So it turns out to be the product of 

(1) a rotation from the fixed frame to an intermediate frame (i/, y/, h) with the angular velocity 

a>i=^Nx; (15) 

(2) a rotation from the intermediate frame (i/, y/, h) to the orbital frame with the angular 
velocity 

G 

<o K = — h. (16) 

r z 

In the absence of perturbations, (o = cjk- This fact led Hansen to call the base (i/, yv, h) 
an ideal frame. The direction xj is called a departure point in the orbital plane. 
Let 6 be the angle from the departure point x/ to the position vector i, so that 

x = xi cos 0+ ji sin 0. (17) 

Then, according to definition (14), the attitude of the ideal frame with respect to the fixed frame is 
described by the differential system 



— xi = coiX xi = — — Nh sin 0, 
at ir 



A 

— ji = (oi X yi = — Nh cos 0, 

dt ir 



(18) 



(19) 



7 

— ti = o)i X h = — N(xi sin 6 — ji cos 0). 

dt G (20) 

The system is exempt of singularities for the class of Keplerian motions that are not linear (i.e., 
for which G = 0). But it is redundant since the unknowns i/, y/, h must satisfy at all times the 
invariant relations 

xi ' xi = ji - yi= h • h = 1, (21) 

xi ' yi= yi ' h = n • xi = 0. (22) 

Nevertheless, Andoyer [12] suggested to integrate numerically the first two equations while 
reserving the third one and the invariant relations to monitor the propagation of round off and 
truncation errors in the course of a numerical integration. 

One way of eliminating the liaisons is to represent the rotation by its Eulerian angles. Let 
h be the longitude of the node of the instantaneous orbital plane in the orbital plane at epoch, / 
the inclination of the instaneous orbital plane over the orbital plane at epoch, and a the argument 
of the departure point reckoned from the node in the instantaneous orbital plane. Note that — cr 
is what classical authors call the longitude of the departure point. The following relations can be 
found in textbooks in kinematics: 

CO/ • ki = h sin / sin cr + / cos o", (23) 



o>/ • ji = h sin / cos a — I sin cr, (24) 

= CO/ - h = h cos / + o\ (25) 

According to (25), once h and / have been determined, the argument cr results from a quadrature; 
in other words, it is determined only to the addition of an arbitrary constant. Therefore, one is 
free to select the departure point at epoch. The obvious choice is the direction x which makes 
the ideal frame at epoch coincide with the orbital frame at epoch. 

Substituting definition (15) in the relations (23), (24) and (25) provides the differential equations 
in the Eulerian angles: 

A = — N sin (0 + or) cosec /, (26) 

1 =-/Vc<>s (0+<r), (27) 



The initial conditions are 



-=-Acos/. (28) 



h = /.» = (To = 0, (29) 



which define precisely the point at which the system fails. 

Various combinations of /? and / have been proposed to eliminate the singularity 1 = from 
the system (26), (27), (28). They derive in principle from the selection 

k\ = sin / cos h, k 2 = sin / sin h 

made by de Sitter [13]. The resulting equations are usually given a form best suited to hand cal- 
culations by reference to a table of logarithms or by means of a desk calculator. 

The instantaneous position of the ideal frame may be viewed as the image of the ideal frame 
at epoch by & finite rotation. 

Rodriguez [14J represents a finite rotation by a vector Q: its direction is the axis of rotation, 

its norm Q is set equal to 2 tan — \, where xis the amplitude of the rotation. 
If co is the angular velocity of the rotation, then 

Q=o> + - 2<a xQ+ |(Q-c»)Q. (30) 

Applied to the rotation vector Qi of the ideal frame, the basic identity (30) yields the differential 
vector equation 



0,-£jv 



x + \xXQ l + -^(Q,-x)Q l 



(31) 



It is equivalent to either the differential system (18), (19), (20) or the scalar differential system (26), 
(27), (28). At epoch the finite rotation is the identity mapping, hence the set of initial conditions 

Q/(ro) = 0. (32) 

5 



Equation (31) may be decomposed in the moving frame of reference. Since 

0/ = T. Qi + <»> x Qi 

ot 

where, as usual, — denotes the time derivative in the moving frame, the vector eq (31) becomes 
ot 



s*-h" 



i-JixQ,+] (Qi-x) Q, 



(33) 



Whether, in the fixed frame or in the moving frame, the vector of finite rotation decomposes in 
the sum 



Qi = Qixo-\-Q 2 yo + Q3n = (?ii/ + (? 2 yi + Q*h, 

and the vector eq (33) gives rise to the scalar system 

(?! = g N I cos - - Q 3 sin 6 + - (Q t cos + # 2 sin 0) Q x 

r r 1 1 

O2 = g AM sin + - (? 3 cos (9 + - «?j cos (9 + (? 2 sin (9) Q 2 



Q> = g n 



-- (Q 2 cos 0-<?i sin 0) 



+ J ((?! cos 6 + Q 2 sin 6) Q 3 1 . 



(34) 
(35) 



(36) 



There results readily the invariant relation 



1 A 



<?l<?2-<?l(?2 + 2<?3 = 



(37) 



to be used to monitor the integration of (31). It is strictly the correspondent of the quadrature (25). 
Equation (31) is not used in VOP, for it can be simplified to the point where its programming 
becomes economical. Starting from Rodriguez identity 



x — Xq = - ()/ X (x + i ), 



one obtains in turn that 



(38) 



x = - (x + Xo) +- 0/ X (x + Xo), 



xXQ,= ±(x + Xo) X Q/ + -|^ (* + *•) —|Q/- (* + «o) Q/, 

1 



0/ ' x = - 0/ • (i + i ) 



(39) 
(40) 

(41) 



On substituting (38), (40) and (41) into (30), one obtains eventually that 

2 0/ = £#(1 + |q?) (x + xo), 

which is one of the equations in Herget's VOP. 



(42) 



The differential equation for the vector of finite rotation is exempt of singularity at 7=0. But, 
should the finite rotation approach a half turn, the amplitude x would come close to 7rand the norm 
Qi would grow beyond any finite bound. Equation (31) fails for half turns of the ideal frame. A sort 
of rectification should remedy this defect. Would the norm of Q/ becomes too large to be significant, 
the last obtained ideal frame should be taken as the fixed frame of reference. However easy in 
principle, this resetting of the epoch imposes various detailed adjustments in the programming 
of the VOP. We favor replacing the vectorial differential eq (39) by one that is definitely free from 
any sort of geometrical or kinematical singularities. 

To this effect we choose to represent the vector of finite rotation Qihy its Eulerian parameters, 
namely 

Xi= (Q/ * x ) sin - x=sin -/-cos -(A — cr), (43) 

X 2 = (0/ * yo) sin - x = sin - I • sin - (A — a) , (44) 

X 3 =(Q/-n) sin -x = cos-/-sin - (A + cr), (45) 

X cos ~x = cos-/ • cos- (A -ho-). (46) 



As we remark that 



Xgd +J<??) = 1, (47) 



Q/-0/=U+^<? 2 ,) (Q/- «/), (48) 



we infer that 

Xo = ^X (Q/ • co), (49) 

-f \ Qi = k (»/ + \wXQj). (50) 

at I 

There remains now to replace Xo 0/ by its development 

XoQ; = 2(Xii/+ \ 2 yi + X 3 h) 

in (50) and it will turn out that the vector eq (33) is equivalent to the scalar system. 

2 \, = /V (Xo u-\ 3 v), (51) 

2 k 2 = N (Xo v+ Kmi), (52) 

2 X 3 = N (X, t;-X 2 "), (53) 

2 ko = N (-Xi u-X 2 v). (54) 



We have just introduced the notations 

" = £ C0S *' (55) 

^i™ ' ( 56) 

to make the equations more concise and to underscore the places where the programming can 
be made efficient. The initial conditions are 

\o=l, X, = \ 2 = X 3 = 0. 

The equations are free from singularities, either geometric or kinematic. Moreover, the Eulerian 
parameters satisfy the inequalities 

-1 ^ \i i^ + l, ^ i^ 3, 

a circumstance that is especially favorable to a numerical integration. The counterpart of the 
quadrature (25) and of the invariant relation (37) is the liaison 

Ai A 2 A. 2 A.] ~~r * Ao A3 Ao A3 — U 

which we readily infer from eq (51) — (54). In addition we have the identity 

\|+\? + \! + A!= 1 

that is inherent to any quadruple of Eulerian parameters. 

Eulerian parameters are somewhat unusual in celestial mechanics. They were considered 
first by Musen to describe the rotation of the apsidal frame [15]; the suggestion was then applied to 
the ideal frame [16] and followed by Broucke for the orbital frame [17]. 

3. The Osculating Ellipse in the Orbital Plane 

Referred to the instantaneous ideal frame, the original equations of motion (1) and (2) become 

x — a>/ X x = — x = X, (57) 

at 

X-o>/XX=-X = -4x + /?i+ry. (5 8 ) 

dt H 

Likewise the equation for the angular momentum converts into 

^G=rTh, (59) 

at 

and the rotation of the orbital frame with respect to the instantaneous ideal frame is described by 
the equations 

— x = 6> K X x = — y, (60) 

dt r l 

8 



h'-***'—?*' (6D 

— ri = (o K X n = 0. (62) 

a? 

Whereas the concept of angular momentum is the root of our considerations about the rota- 
tion of the orbital plane, the Laplacian vector 

A = XxG — fix (63) 

is the cornerstone in our description of the osculating ellipse about the center of attraction in the 
orbital plane. The dominant role in the dynamics of the problem is played by the effective perturbing 
force 

K= F-- ix (ixF), (64) 

P 

where 

P V (65) 

is the semiparameter of the ellipse. In the sequel the following identities 

K- x = F •*, 
K-y=(l + r -)F.y, 

K - A = 1 + - ) F • h 



will be used repeatedly without being referred to specifically. 

Cross multiplication of (63) by h yields the equation of the hodograph 



X = ^UxA + gj 



and brings us to consider the vector 



B=-nXA (67) 

in the direction of the minor axis. 

On introducing definition (67) into the equation of the hodograph (66) and then differentiating 
with respect to time in the instantaneous ideal frame as per eq (58) and (61), one finds eventually 
that 

-B = K-(K- h) h= (hxK) x h. (68) 

at 

Let b be the direction of B; also let /ne denote the norm of the Laplacian vector. Then (67) takes the 
form 

G b (69) 

9 



and, on using (68), one finds readily that 

dt G ' (70) 

which is the most concise expression we could find for the time derivative of the eccentricity. 

If a denotes the direction of A, then h = a X b and the triple (d, b, h) constitutes that ortho- 
normal base that we proposed in the introduction to call the apsidal frame. With the help of (68), 
(69), and (70), one checks readily that 

Q 

(*>a = (K • a) h (71) 

fie 

is the angular velocity of the apsidal frame with respect to the instantaneous ideal frame, so that 
the rotation of the apsidal frame is described by the equations 

A C 

— a = co., X d = (K • a) b , (72) 

dt jxe 

r) » C 

— b=o) A x b= — (K-d)a, (73) 
dt fxe 

— h= a) \ X h = . 

dt 

The argument of perigee g is defined here as the angle from the departure point xi to the semi- 
major axis a. In order to circumvent the indetermination due to small eccentricities, we consider 
the functions 



/ne 
C = -cosg, (74) 

S = ^sing. (75) 

The equations for C and S result from differentiating the sum 

Lr 

in the ideal frame to obtain that 



dii+ ^ ,= di{f)- d+t fk d=Kx ^ 



thus leading to the simple (and probably original) equations 

C = Kyj, (76) 

S=-Kxi. (77) 

The motion of the osculating ellipse in the orbital plane with respect to the instantaneous ideal 
frame is completely determined by the system consisting of eq (9), (76) and (77). 

10 



Special care must be taken while computing a position on the osculating ellipse in order not 
to bring in singularities (e.g., divisions by e) or indeterminations (e.g., appearances of sin g and 
cos g in other combinations than C and S) caused by small eccentricities. We postpone until the 
next section dealing with the problem of timing the motion along the osculating ellipse, and so 
assume that the mean longitude 

F = l + g (78) 

reckoned from the departure point is known. Once C and S have been obtained, we revert to the 
more commonly used elements 

X = ^C = ecosg, (?9) 

We also consider the eccentric anomaly 

<l) = E + g (81) 

reckoned from the departure point. The classical Kepler equation 

I = E — e sin E 
may be extended into an implicit equation relating </) to F: 

F = <f> - (X sin </> - Y cos (/>). (82) 

Likewise we transform the classical identities related to the anomalies reckoned from the line 
of apsides d 

r = a (I — e cos E), 
r cos /= a (cos E — e) , 
r sin /= a r) sin E, 
in which 

1) = (i-e 2 ) 1 / 2 , 
into formulas involving 6 and </> reckoned from the point of departure: 



- = 1 — X cos (/> + Y sin </), 
a 

<f>-F 



(83) 



^cos0 = cos4>-X + f— L F, (84) 

a 1 + t? 

11 



with 



- sin 6 = sin <b — Y : — ; X, 

a 1 + 7} 



v = (!_^2_p)i/2 



W 



(85) 

(86) 
(87) 



4. The Timing on the Osculating Ellipse 

Most authors borrow the differential equation for F from Brown [18]. His derivation rests on 
the use of a certain differential operator which he denotes 5; the present exposition shows that it 
is not needed. Also Brown's construction dates from a time when astronomers were not yet aware 
of vector calculus. 

We start with writing the time derivative of the mean anomaly as follows: 



l = E r (e sin E) ; 

dt 



(88) 



each term in the right-hand member will be calculated separately. 

Introducing the equation of the hodograph (66) into (11) results in the expression 

r = — sin / 



(89) 



Now, by definition of the eccentric anomaly, 

r sin/= a r) sin E\ 

hence, formula (89) takes the form 

r r = L e sin E 
provided that 

L= (/xa) 1 / 2 



(90) 
(91) 



denotes Delaunay's action variable. Actually (90) yields the time derivative of e sin E as the 
combination 



d _ 1 

— e sin L = — 
dt L 



l (rr) - rr i logL ]- 



We evaluate 



^rr=^~ (x -X) = (X-X) + ( x ^X 
dt dt \ dt 



= (X-X)-^+ ( x . *•)=-■£ + ■£+ (x-F) 



and obtain eventually that 
d 



. e sin E = — n + — (K • x) -f 
dt L 



a r r d , 



(92) 
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In order to compute E, we start from the identity 

Ed, E 

^—F = ~L log tan - 
sin E at 2 

which, in view of the classical formula 



transposes into the identity 



E l-e f 

tan 2 = A IT7 tan 2' 



L^f 



sin E r\ z sin/ ' 

Reduction of the first term in the latter sum goes as follows 



e e 2 + ri- , / d , 

— e = e \—r log e 



T) 2 Tj 1 \dt 



i logr >) 



(93) 



On account of (70), the reduction may proceed further to 

^ = - (K-5)+ e ( T log L. (94) 

i)- i± at 

As for the second term in the right-hand member of (93), we observe that 

f=6-g, (95) 

. G 

= -ii (96) 

r- 

G v m 

g= K- a = €0 A 'n. (97) 

fie 

Equation (96) is a consequence of the hodograph eq (66). Indeed from taking the time derivative of 
x in the instantaneous frame of reference, we obtain immediately that 

so that 

r = X • y=£ (1 + e cos f)=yr = — . 
G Gr r 

As, for eq (97), it derives readily from differentiating the definition (74) on account of (76). 
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By combining expressions (94) — (97), we expand (93) into the sum 
Before we combine (92) and (98) we observe that 



/HI 

E= K*(d — besinf) + 

\xar\e 



a — be sm f = — a — ex. 

r 

Eventually we reach a rather concise equation for the mean anomaly: 

/=ri+^(K-d)-2T(K-*). (99) 

jxe L 

As expected, it is singular at e = 0. The critical point is eliminated by adding (97) and (99), which 
produce the equation for the mean longitude 

F = n+ I (K-d)-2J(K-x). (100) 

Eventually to avoid having to determine the vector d, we decompose the right-hand member of 
(100) in the instantaneous ideal of frame, which produces an equation exempt of any problems 
for small eccentricities: 

F = n+ /n P - (CS-CS) + 2r)(uS-vC). (101) 



5. Conclusions 

Completing the project opened by Herrick, developed by Musen and revised by Herget, 
table 1 summarizes the final version of a VOP free of geometrical and kinematical singularities 
and indeterminacies, with the right-hand members prepared for as efficient a coding as possible. 

The elements (Xo, Xi, X2, X3, G, C, S, F) are termed ideal because they refer to an ideal frame 
of reference. 

The derivation presented in this communication departs from the general methods of the theory 
of perturbation based on the Lagrange brackets for the conventional elliptical elements [19]. Rather 
it is based on the method of variation applied immediately to the vector integrals (angular 
momentum and Laplacian vector) of the problem of two bodies and the well-known formulas for 
differentiation in moving frames of reference. We concur with Lure [20] and Danby [21] in thinking 
that this is the most simple, informative and economical of all extant derivations. 

Table 1. Variations of parameters 

Attitude of the Ellipse in the 

ideal frame orbital plane 

2\ l = N(\ou-\ 3 v) G=rT 

2k 2 = N(k t) v-\-\3u) C = Ky, 

2ka = N(k l v-k 2 u) -S = K-x, 
-2ko = N(kiU+k 2 v) 
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Timing 



F=n+ — -^ {CS - CS ) + 2i } (uS - vC) 

/J. ( 1 + 7/ ) 

Invariant relations 



\J + \f + A^ + Xr; = 1 
\o\3 — XoAa + XiA* — A 2 Xi = 
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